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Abstract 

We study aspects concerning numerical simulations of Lattice QCD with 
two flavors of dynamical Ginsparg-Wilson quarks with degenerate masses. A 
Hybrid Monte Carlo algorithm is described and the formula for the fermionic 
force is derived for two specific implementations. The implementation with 
optimal rational approximation method is favored both in CPU time and 
memory consumption. 
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1 Introduction 



Recently, in a series of publications |[ ^, ^, ||, ||, it has become clear that, if one would 
modify the chiral transformation laws away from their canonical form in the continuum, chiral 
symmetries can be preserved on the lattice without the problems of fermion doubling. The 
lattice Dirac operator D for these fermions satisfies the Ginsparg- Wilson Q relation^], 

^5D + D^5D = aDj^D , (1) 

where a is the lattice spacing. As a consequence of the Ginsparg- Wilson relation (^, it is easy 
to show that the fermion action, 

5f = , (2) 

x,y 

is invariant under lattice chiral transformations, and chiral symmetry will protect the quark 
masses away from the additive renornializations. 

The chiral properties of the Ginsparg- Wilson fermions is a direct result of the Ginsparg- 
Wilson relation. In particular, several types of fermion actions could be written down which 
all fulfill the condition. For dcfinitcncss, one particular choice ^ is adopted in this paper, 
namely: 

aD = 1~Hw{hI^Hw)-^'^ , 



Hw - l + *--a^^[7M(V^ + V;)-aV;V^] , 

= (-3 + s)(5,,j, + i^ [(l-7,0C/p(a:)W!/ + (l + 7M)t^^(a;-A*)'5..-M,y] (3) 

where s is a parameter satisfying \s\ < 1. The lattice covariant derivatives V* and are 
defined as usual, 

V;V = ^P{x)-Ul{x-^,)i;{x) . (4) 

Due to their decent chiral properties, it is quite tempting to investigate the possibility of 
performing Monte Carlo simulations using Ginsparg- Wilson fermions, replacing the conventional 
Wilson fermions. Albeit their seemingly non-local appearance, the fermion matrix is in fact 
local (with exponentially decaying tails), as long as the parameter s in the matrix is chosen 
appropriately ^ . The locality property of the fermion matrix will enable us to use iterative 
methods in Krylov space whenever inversion of the matrix becomes necessary. 

When performing the inversion of the fermion matrix (H) , one would encounter the problem of 
yet another matrix inversion of a fractional power. Recently, proposals have been put forward 

^They are therefore named Ginsparg- Wilson fermions. 
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^ which make such inversions possible. Some numerical calculations in quenched QCD 
[ pi] for Ginsparg- Wilson fermions already indicate that these fermions indeed have anticipated 
chiral properties. However, it is also realized that quenched calculations with Ginsparg- Wilson 
fermions are more costly than with conventional Wilson fermions primarily due to fractional 
inversion of the matrix H^^Hy/- Two specific types of methods will be discussed in this paper. 
One is the method proposed by Bunk p^ , which will be called fractional inverter method, or 
FIM. The other method is the optimal rational approximation method, or ORAM, proposed 
in Ref . . It was reported in Ref. that ORAM converges faster than FIM for a desired 
accuracy and a given condition number of the matrix H^Hw- Now, each multiplication with 
the fcrmion matrix D for Ginsparg- Wilson fermions is equivalent to 2N + 1 multiplications 
with matrix i/^ or Hw, where N is some integer. With FIM, N is the number of highest 
order Legendre polynomial kept in the iteration procedure. With ORAM, TV = Ncg is the 
number of conjugate gradient |^ iterations needed to perform the multi-shift matrix inversions. 
This number is determined by the condition number of the matrix H^^^Hw and the accuracy 
desired. Therefore, the calculations with Ginsparg- Wilson fermions is at least more costly than 
conventional Wilson fermions by a factor of 2N + 1. 

Although it is already quite costly in the quenched case, it remains an tempting problem 
to simulate dynamical Ginsparg- Wilson fermions. No algorithms including dynamical fermions 
have been tested on these newly proposed fermions. In this paper, it is shown that a Hybrid 
Monte Carlo algorithm would do the job, however, just as in the quenched case, the simulation 
is more costly than simulating dynamical Wilson fermions. Also, in the calculation of the 
fermionic force to gauge links, using two different methods for the fractional inversion results in 
very different memory and CPU time consumptions. For FIM, it seems that 0{N) psudofermion 
fields have to be stored in order to make the simulation tractable. For ORAM, only a moderate 
number of psudofermion fields have to be stored, and the number of matrix multiplications 
also increases slower than in FIM. In this paper, we will concentrate on a Hybrid Monte Carlo 
algorithm for simulation of dynamical Ginsparg- Wilson fermions. The fermionic force for the 
gauge links, which is the crucial part for the dynamical fermion simulation, will be calculated 
for both ORAM and FIM. General properties of the two methods are compared. Test runs on 
small lattices with gauge group SU{3) are now being investigated and detailed results will 
be reported later. 

This paper is organized in the following manner. In Section 2, the Hybrid Monte Carlo 
algorithm suitable for simulating dynamical Ginsparg- Wilson fermions are described and the 
formula for the force of the gauge link is derived, in both ORAM and FIM. These two methods 
are compared in the dynamical simulation in terms of CPU time consumption and memory 
consumption. Possible improvement methods arc also addressed. Some concluding remarks are 

^ By the phrase "conjugate gradient", we mean all possible iterative algorithms in Krylov space: conjugate 
gradient, minimal residue, Bi-conjugate gradient, etc. 
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in Section 3. 



2 The Hybrid Monte Carlo algorithm 

The basic formalism of Hybrid Monte Carlo algorithm jljl remains the same as in the conven- 
tional Wilson case . Only the fermionic force has to be re-derived for the Ginsparg- Wilson 
case, which will be dealt with below in detail. The effective action with the psudofermion 
contribution now reads: 

S,ff = Sg[U^{x)]+cf>^Q-^cl>, (5) 

where the fermion matrix Q = ^^{D + m) is hermitian and (f>{x) being the psudofermion field 
generated at the beginning of a Hybrid Monte Carlo trajectory from Gaussian noise. We have 
also assumed that two flavors of quarks have degenerate masses. At each molecular dynamics 
step in a Hybrid Monte Carlo trajectory, one has to find solution vectors Xi = Q~^(f> and 
X2 = Q^^4> — QXi from an iterative algorithm (conjugate gradient, for example) in Krylov 
space. Then, the total force with respect to the gauge fields can be found by investigating the 
variation of the action under infinitesimal changes of the gauge fields: 

6Sf = S[(l)''Q-^(l)]^Tr[F^{x)6Uf,{x) +h.c.] . (6) 

The gauge staple Vfj_{x) comes solely from the pure gauge part of the action and could be 
obtained with little cost. The fermionic forces Ffj_{x), however, is much more costly. Once the 
fermionic forces are obtained, the standard Hybrid Monte Carlo updating procedure can be 
carried on just as in the conventional Wilson case. 

To derive the formula for the fermionic force, we take the variation of the fermionic part of 
the action and get, 

SSf = Xl{-SQ)X2 + xl{-SQ)Xi (7) 
The variation of the matrix Q contains two parts, one being simple, namely 

-SiQ = j,{SHw)iHlHw)-'/^ , (8) 
the other being quite complicated, i.e. 

-52Q^l5Hw5{H\vHwy^'^ , (9) 

which depends on the detailed implementation of the fractional inversion of the matrix. We 
now proceed to discuss the fermionic forces in ORAM and FIM respectively. 
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2.1 Fermionic force in ORAM 

We first present the force in ORAM wliicli is more straightforward. We recall that, this approx- 
imation amounts to approximate the function z(z^)^^/^ in the interval [0, 1] by a ratio of two 
polynomials: 

.(.r"^=.(c.+|;^) , (10) 

It is an approximation similar to the Fade approximation used in Rcf . [p6| . For details about 
this approximation and how to get coefficients Ck , consult [ p] and references therein. Applying 
this method to the hermitian matrix j^Hyy, we immediately obtain the following expression for 
the variation of the fermionic action: 

N 



SSf = coTriX2<i^Xh5SHw) + Y.CkTriC2,k®Xlj5SHw) 

k=l 
N 

+ CoTr{Xi ® xl-f^SHw) + ^ CkTr{Ci,k ® xIj^SHw) 

k=l 

N 

k=l 
N 

- CkTr (CiM ® ^^k^lvSHw + ^i.k ® 4^kl5SHv 



k=l 
N 

m /a „ rrt err , ^ „ c r, 

'TV 

fe = l 

Q,k = 7 „ ^2 I ^» ' ii,k ^ l^Hw-, „ .2 , ■ (11) 

KlbHw) +qk hbHw) + Qk 

In the above formula, "Tr" stands for taking trace in both Dirac and color space and a sum- 
mation over the whole lattice points. The symbol (E> stands for direct product of two vectors in 
color space. Therefore, the fermionic force is obtained as 

F^{x) = ^■trD^rac[X2{x + fi) (E) x|(a;)75(l - 7^) + Xi{x + /i) (g) A:2^(a;)75(l - 7^)] 

N 

Ck , 



+ ^ trmracY[C'2,k{x + fi) <Sl X|(x)75(l - 7^)] 
k=l 
N 

+ ^ i?-£)i™cy [Ci,fc(a; + A*) ® Xl{x)j5{l - 7^)] 
fc=i 

N 

- ^troirac-^lGA^ + A^) «> [Hw^i,k]''{x){l - 7^)] 

N 

- ^ trDtrac — [S.2,k{x + A*) «> (a;)75 ( 1 - 7p)] 
k=l 
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JV 

Cfe 



- y^Jrpirac^iCl.kix + ^) (g) [Hw£.2,k]Hx){l - 7^)] 



k=l 
N 

- ^trDzracYi^iA^ + ® ^iki^h^i^ - . (12) 

fc=l 



where the trace trjjirac is only taken withm the Dhac space. 

2.2 Fermionic force in FIM 

In FIM, we would like to solve for ^ satisfying: 

Af 1/2^ = X , given X. (13) 
where the matrix M = H^^Hw by setting: 

M = c(l + <2 _ 2tA) , (14) 

with the parameters c and t chosen in such a way that all eigenvalues of the matrix A lie within 
[—1,1]. To be more specific, we choose, 

o.<4±i£, (1.) 

where Xmm i^max) is the lowest (highest) eigenvalue of the matrix HlyHw and k = A 

max I ^min 

is the condition number. Then, the solution to eq. ( p^ ) may be written as: 

oo oo 

e = c-i/2^t'=P,(^).X = ^Sfc , (16) 

fe=0 fe=0 

where Pk{z) are Legendre polynomials. Therefore, an approximant for the solution at the n-th 
level is 

n 

in=Y,Sk . (17) 

fc=0 

The shifts ^, Sk, defined as 

Sk = c~^'h''Pk{A)-X , (18) 

satisfy the following recursion relations: 

s_i ,so = c^^/^X , 

Sk = {2^1/k)tAsk-i-{l-l/k)t^Sk-2 ■ (19) 

•^A extra factor has been included in the definition of as compared with Ref. |12| 
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For the case of Legendre polynomials, it is claimed that the following bound for the residue is 
obtained 111: 



le - urn = \Rn{A)\ < = (^^)"+^ , (20) 



which asserts the exponential convergence of the iteration. 

For the vectors 6{Hl^Hw)~^^'^Xi, similar strategy could be apphed, 

(5(Af-i/2)?7 = c-i/2^t"JP„(A)ry , given rj, (21) 

71=0 

where rj represents either Xi or X2 and Ncut is the highest order of Legendre polynomials kept 
in the approximation. In an analogous manner, we define, 

6k = c-'^h\6PkiA))-r^ , (22) 

which satisfy the following recursion relation: 

{k + 1)4+1 + kt^k-i = (2fc + mSAsk + A5k) . 

5_i =0 ,4 = 0, (23) 

Using this relation, 4 could be expressed as: 

fc-i 

5k^Y.*^kiA)SAsi . (24) 

1=0 

The coefficients L''f,{A) are polynomials in A with degree {k — I — 1) and can be expressed as 
Legendre polynomials, 

4(A)^f^«"-n-,-(f^A, . (25) 

After rearranging the double summation and some trivial algebra, we get the following formula 
for the variation of the fermionic action: 

6Sf = Tr{Zi (g) xIj55Hw) + Tr{Z2 ® xIj^SHw) 
1=0 



2c ^ 

1=0 



21 - 


f 1 


U 


- 1 


21 - 


f 1 


U 


- 1 


21 - 


f 1 



/=0 
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1 0/^1 
1=0 

= [HlHwr'^^Xi , Xi,i=t\Pi{A))Xi , 



= E <f =t'"P^((f^)75^75)if^X, . (26) 

rn=i) 

Therefore, the following expression for the fermionic force is obtained: 

Fi,{x) = trDirac{Zi{x + n) Xl{x)'y5{l - j^,) + {Z2{x + n) Xl{x)'y5{l - 'y^,)) 

/V 

1 J ^ 21 + 1 

~ o~ E( ; I 1 )^^J'i'-ac[(-gtya:2,;(a: + /x) (g) Ti,,(a;))75(l + 7^)] 
1 ^ 2Z + 1 

- 2-E(T:rr^*''^'™'=[^^2''^^ + '^) ® [-ffjyri,i]^(a;))75(l -7/.)] 
1 ^ 2? + l 

- ^ E'^Trr)*''^'™'=[('^^^i''(^ + ® 72,;(2^))75(l + 7/^)] 

/V 

1 2Z + 1 

- ^E(TVY)*^ci™c[K;(a; + M)®[ifi^r2,i]t(a;))75(l-7M)] • (27) 

Since the CPU cost of a simulation program with dynamical fermions is dominated by the 
fermion matrix times vector operations. It becomes clear that the above formula for the 
fermionic force is not very useful from a practical point of view. The most CPU consum- 
ing part is the calculation of the vectors Ti^i{x), for all values of I, each requiring an iterative 
procedure, i.e. the calculation of the quantities S^^K This implies that, in order to calculate the 
fermionic force, the multiplication of the matrix hI^^Hw has to be performed 0(N^^^) times, 
where Ncut can become large. This would make the calculation of the fermionic force too costly. 

However, there is a way around this difficulty. The price to pay will be some extra storage. 
Note that Legendre polynomials Pm{P{l)z), (3{l) being the constant {21 + + 1), could 
be expressed as a linear combination of Legendre polynomials of lower or equal degrees with 
argument changed to z, i.e. 

m 

Pm{Pz) = ^a,n,MPj{^) ■ (28) 

j=0 

With this, we could express the quantities Tj^; in the following way: 

N-l-l 

Ti,l = E fm{l,t)t'^Pm{l5Aj5)HwX, , 
m— 
N~l-1 

fm{t,i) = Yl t'-"''^j,mm)) ■ (29) 

j=m 
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The functions fm(t, I) are just c-numbers and can be calculated at the beginning of the simula- 
tion. Therefore, after the vectors Xi are obtained, one can calculate the vectors Pmi'jsAj^) H Xi 
once for all values of m and store the resulting vectors. Thus, Ti^i could be obtained easily 
without further iteration of matrix multiplications. The coefficients (Tm,j(l) satisfy the following 
recursion relation: 

where the subscripts m j should satisfy < j < m and a zero value is understood whenever 
an out-of-range subscript is encountered. Together with cro,o = Ij the above recursion relation 
completely determines all coefficients am,j{P{l)) and therefore the function fm{t,l)- Now the 
calculation of the quantity Tij only requires a linear combination of vectors which costs little 
CPU time. 



2.3 comparison of the two methods 

We now compare the CPU time consumption and memory consumptions of the two methods 
discussed so far for the simulation of dynamical Ginsparg- Wilson fermions. As is well known, 
the CPU time consumption is proportional to the number of operations of the matrix Hw 
on vectors. For each molecular dynamics step in the Hybrid Monte Carlo, ORAM requires 
2-^VcG(2-/Vcg + 1) number of matrix multiplications to obtain the solution vector Xi and 4A^cg + 
4:Nr more matrix multiplications to obtain the fermionic force. Here, Ncg is the number 
of conjugate gradient iterations needed to obtain the solution Xi and Ncg is the number of 
conjugate gradient iterations needed to obtain the vector {H^Hyy) + qmin)~^ Xi for the smallest 
shift qmin- Parameter Nr is the order of the polynomials in the optimal rational approximation. 
ORAM also requires to store Nr psudofermion field vectors. As a comparison, FIM requires 
2NcG{'^^cut + 1) number of matrix multiplications to obtain the solution vector Xi and 12Ncut 
more matrix multiplications to obtain the fermionic force, where N^ut is the highest order of 
Legendre polynomials kept in the series expansion. FIM also needs to store 2Ncut psudofermion 
field vectors. 

Concerning the CPU time consumption, both method are more costly than dynamical sim- 
ulations with Wilson fermions by a factor of 2N + 1, where N = Ncg for ORAM and N = Ncut 
for FIM. From the theoretical upper bound of the error, the two methods behave in a similar 
manner, Ncg ~ Ncut- Practically, however, according to the experience in [ pT[ , Ncg is usually 
less than Ncut because the theoretical bound is saturated for FIM while it is not for ORAM. 
Therefore, ORAM is more favorable compared with FIM when doing simulations with dynam- 
ical Ginsparg- Wilson fermions. Concerning the memory consumption, N^. is usually much less 
than Ncut which would again put ORAM in a more favorable place. 

It is clear from the above discussion that, if one would like to accelerate the simulation with 
dynamical Ginsparg- Wilson fermions, one has to find ways to decrease Ncg in ORAM or Ncut 
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in FIM. These two parameters are mainly determined by the condition number of the matrix 
h'^Hw- Any preconditioning methods that would decrease the condition number of the matrix 
(while still maintaining the shifted nature of the matrix in ORAM) will bring an improvement 
to the simulation of dynamical Ginsparg- Wilson fermions. It should be pointed out that other 
improvements, for example using better integration schemes, would apply to both methods. 
Test runs on small lattices are now under investigation [ p^ where these algorithmic issues will 
be further studied. 

3 Conclusions 

In this paper, possibilities of simulating dynamical Ginsparg- Wilson fermions are discussed. The 
formula for the fermionic force is derived for two specific implementations of the algorithm, the 
optimal rational approximation method (ORAM) and fractional inverter method (FIM). It turns 
out that, in both methods, simulating dynamical Ginsparg- Wilson fermions are more costly than 
simulating dynamical Wilson fermions. The extra CPU time consumption mainly comes from 
the the fractional inversion of the matrix. In quenched simulations, it has been reported [ p] 
that ORAM performs better than FIM. In dynamical simulations, this conclusion remains, 
both on CPU time consumption and memory consumption. It should be emphasized that, 
though being more costly, the advantage of simulating dynamical Ginsparg- Wilson fermions 
over dynamical Wilson fermions or quenched Ginsparg- Wilson fermions would be a much better 
behavior towards the chiral limit. The feasibility of such simulation has been demonstrated in 
this paper using a Hybrid Monte Carlo algorithm. 
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